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Abstract 

The effects of viscosity and small-scale atomic-level mixing on plasmas in inertial confinement fusion (ICF) 
currently represent challenges in ICF research. Many current ICF hydrodynamic codes ignore the effects 
of viscosity though recent research indicates viscosity and mixing by classical transport processes may 
have a substantial impact on implosion dynamics. We have implemented a Lagrange hydrodynamic code 
in one-dimensional spherical geometry with plasma viscosity and mass transport and including a three 
temperature model for ions, electrons, and radiation treated in a gray radiation diffusion approximation. 
The code is used to study ICF implosion differences with and without plasma viscosity and to determine 
the impacts of viscosity on temperature histories and neutron yield. It was found that plasma viscosity 
has substantial impacts on ICF shock dynamics characterized by shock burn timing, maximum burn 
temperatures, convergence ratio, and time history of neutron production rates. Plasma viscosity reduces 
the need for artificial viscosity to maintain numerical stability in the Lagrange formulation and also 
modifies the flux-limiting needed for electron thermal conduction. 


1 Introduction 


Direct-drive inertial confinement fusion (ICF) refers to laser heating of a spherical shell that contains fusion 
fuel, resulting in compression and fusion 03 Perfect spherical compression has many shortcomings resulting 
from drive asymmetries, laser-plasma instabilities, and fluid instabilities including Rayleigh-Taylor (RT), 
Richtmyer-Meshkov (RM) and Kelvin-Helmholtz (KH). Target capsule non-uniformities as well as beam- 
to-beam imbalances become exaggerated during implosion convergence and perfect spherical compression is 
compromised and degrades performance. The implosion dynamics can also be affected by the shell material 
mixing into the fuel, and the mix can in turn be influenced by plasma viscosity. Plasma effects on burn physics 
and implosions have been the focus of many studies, including work on the plasma mix layer structure [3j, 
barodiffusion [4], flux limiters §0. as well as laser absorption mechanisms [8| and asymmetries [9 mostly 
using hydrodynamics codes in comparison to experimental data. 

There has been increased recent interest in plasma kinetics as a possible mechanism to explain aspects 
of ICF implosion degradation 10 m- Kinetic theory describes plasma transport in the small Knudsen 


number limit resulting from particle collisions and attempts to determine consistent transport coefficients 
multiplying the known gradient quantities to determine species diffusivity, viscosity and thermal conduction 
for the mass, momentum, and energy conservation equations 12-14 . In ICF simulations, energy transport 
typically includes the classical species temperature coupling and thermal conduction [i.e., Spitzer 


15 


Conduction is dominated by the electron contribution and is often flux limited as a mechanism to improve 
simulation results in comparison to experimental data |6}|7]. Species plasma viscosity is usually assumed 
to be negligible [l 16 and may be overwhelmed by artificial viscosity 17 used in many Lagrangian-based 
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hydrodynamic simulation codes for numerical stability and to capture the shock discontinuities. Particle 
diffusion is often assumed to be negligible while species mixing is modeled as the result of classical fluid 
instabilities (RT, RM, KH) [I] or from turbulent processes 18 19 . Recent analysis 20 suggests that the 
mix in ICF implosions involves comparable contributions from chunk mix (e.g., subgrid scale fluid structures) 
and atomic mix, suggesting that plasma transport, in viscosity and in species diffusion, may have significant 
contributions to mix. 

While previous studies examine various physics issues that may affect implosion dynamics, few specifically 
address the effect of plasma ion viscosity on ICF implosions. Early on, Yabe and Tanaka 21 reported that 
plasma viscosity is likely to be important in ICF implosions, and following up on that, simulations confirmed 
that kinetic effects may increase the viscosity at high Mach number shocks 22 . In a more recent study, 23 


post-processing of ICF simulations was used to infer that viscous stress is small compared to the plasma 
kinetic pressures. In another study, post-processing ’Omega-scale’ ICF simulations, 24 , the time integrated 


viscous diffusive lengths were estimated to reach the dimensions of the compressed fuel during the maximum 
compression. The differences and the uncertainties in transport estimates made by post-processing Eulerian 
simulations highlight the sensitivity of the results to the dynamic trajectories of the plasma densities and 
temperatures while also implying a need for self-consistent in-line plasma transport calculations. 

More recently, the role of plasma viscosity in smoothing small scale fluctuations in 2D ICF simulations 
has been reported 25 , where a distinct smoothing of small scale structures is evident in the ICF fuel region 
when viscosity is included. Plasma viscosity and diffusion have been shown in 2D simulations to play a 
significant role in attenuating RT and KH interfacial instabilities for mode scale lengths and conditions 
relevant in ICF 26 . These studies suggest that plasma dissipative processes are important at micron scale 


lengths for ICF relevant conditions. This in turn suggests that high resolution 3D ICF simulations 18 19 


that now solve the Euler equations on sub-micron scale grids should include the plasma dissipative processes 
in the simulations for accuracy at the smallest scales. 

In a study examining the physical plasma viscosity 27 and comparing to Von Neumann artificial viscosity 


17 , it was found that the plasma viscosity in 2D geometries and conditions relevant to ICF is likely to be 


important as viscous dissipation can increase the fuel entropy and degrade implosion performance. Similarly, 
our present study implements plasma viscosity in addition to artificial viscosity. We found that the computed 
results remain numerically stable if artificial viscosity is turned off after the first shock converges on center. 
The solutions, with artificial viscosity set to zero after the first shock convergence, lead to greater fluctuations 
in density and in temperatures. Assuming these fluctuations reduce the net radial implosion velocities, this 
implies one might expect a decrease in predicted performance using Lagrange codes with no artificial viscosity. 

In the present study, we explore the dynamic effects of viscosity on ICF implosions through the use of 
a one dimensional, three temperature, Lagrangian hydrodynamics model that also includes a treatment of 
fuel-plastic mass mixing by plasma transport. The study focuses on ’Omega-facility’ scale ICF experiments 
with simplified capsules of a CH shell and a deuterium fuel region. It is shown that the simulation results 
are generally consistent with experiments and previous simulations. We then examine A-B comparisons 
of simulations both with and without plasma viscosity in the hydrodynamic equations. We demonstrate a 
significant influence on timing in ID implosions, particularly at first shock convergence and burn conditions, 
while yield and time-integrated burn-weighted temperature show less sensitivity to viscosity. 


2 Theory 

2.1 Hydrodynamics 

The plasma is modeled with one-dimensional hydrodynamics equations for the mixture averaged quantities, 
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where 7 = 5/3 for the plasma as a monatomic ideal gas, and n is the number density for particles including 
ions and electrons. Equations ([Tj) to ([3| express the conservation of density and momentum with the evolution 
of temperature, respectively in the Lagrange frame. Note that the latter two equations have diffusive-like 
terms arising from the plasma transport, a viscous stress tensor in Eq. (J2| and the heat conduction in Eq. 
©■ The last term in Eq. © represents the viscous heat dissipation. Enthalpy carried by the species mass 
flux contributes to the energy within the narrow fuel-plastic mix layer formed by mass diffusion, however, 
this heating term is ignored in this study focusing on the viscous effects. 

In the fuel, the transport coefficients are appropriate to a simple singly charged ion plasma 28 , r/ 0 = 
0.96 nkTri and K e = 3.2nkTr e with the collision times as given in 28 29 . The viscosity is modified in the 


higher z capsule material and in the binary mixing layer between fuel and the CH capsule, consistently with 
the low z - high z kinetic formulation in 12 . The approximate viscosity in the mix layer, with subscript 1 


representing the fuel ions, and 2 the average ion in the CH, is evaluated as 


_ rijkTj ~ rijkTi ~ rijkTj ~ mkT f /2 

710 Vi C„(v 1 , 1 + 17 , 2 ) C v (ni+n 2 z%) 

where z\ = 1 and rij = ni+n^- As a further approximation, the coulomb logarithm is assumed to be fixed for 
each pairwise collision rate, L ss 5, and is folded into the coefficient, C„. With the collision rates proportional 
to fcTj 3 ^ 2 , the plasma viscosity scales as kT ^ 2 , and the kinematic viscosity scales approximately as kT^ 2 /rii, 
highlighting the sensitivity of the dynamically evolving viscosity to the local and instantaneous plasma ion 
temperature and density. 

A three temperature model described in a following section provides a more detailed description of the 
plasma. As such, Eqnj3]is modified to include separate temperatures for the three components in the plasma: 
ions, electrons, and radiation. Thermal conductivities for ions or electrons in the 3T model scale as T 5 / 2 
similarly to the ion viscosity. 


2.2 Species Mass Mixing 

In addition to the single fluid hydrodynamics, the code also tracks species mass composition. The equation 
for the mass fraction, Y, for the light ion species, 1, is 

DY 

-j^r = • Pi w \ (5) 

where w\ = u\ — u is a species drift flux, the species average velocity relative to the mixture mass averaged 
velocity, u. The plasma species mass drift flux in binary mixing due to the concentration gradient and 
barodiffusion terms is written for the light ion species, 1 , as 
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where Y is the mass fraction of the light species, pi a , is the total ion pressure, P 2 is the mass density of the 
heavier species, 172 is a momentum exchange rate between the two species, x is th e number fraction of the 
lighter species, and the charge fraction is Z = n\Z\/n e . The transport coefficient, p = (PiaP 2 )/(pv\ 2 &\ l [x\) 
can be shown to be equivalent to that derived in 12 with the kinetic correction factor, af, [y]. The molar 
concentration gradient is converted to a mass fraction gradient to allow an implicit solution to the diffusive 
component of mass flux, using Vy = (Dy/DF)VT, where D\/DY is a function of only the particle masses. 
The form for the first three terms in Eqn. ©. the molar gradient, the ion barodiffusion, and electron 
barodiffusion, are in general agreement in the literature [TTJ 12j|30 . 
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The thermal gradient forces, the last two terms in Eq. have recently been defined in the literature in 
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The thermal force coefficients here, a^ e ,a^ i can be evaluated in terms of the 


varying forms 

kinetic transport coefficients, an,a.T,oi \ 2l given in 12 , which modify the coefficients for the concentration 


varying across the mix layer. Although the coefficients of the temperature gradient force can appear to be 
large 31 our preliminary analysis suggests that when realistic temperature gradients are accounted for, then 
the forces are small in our test problems. 


2.3 Plasma Pressure 

The plasma pressure utilizes an ideal gas approximation p = nkT , appropriate to multi-species, as a sum 
over ion species assumed to be fully ionized ideal gases and assuming ambipolar electrons. The pressure can 
thus be expressed as a function of the total density, temperature, partial ionization, Zj. , and the light ion 
mass fraction, Y. For a single temperature and binary mix of ions with light ion atomic mass, A\. and heavy 
ion atomic mass, A 2 , the pressure is 


p = nkT = (m + n e )kT = m(l + 21 ) + n 2 ( 1 + z 2 )kT = ^-|j-(l + z i) + ~^(1 + - 2 )^ pkT (7) 
and for the 3T model with T e 7 ^Tj, 

p = UikTi + n e kT e = p ^-^-(fcl) + 2 ifcT e ) + kT t + z 2 kT e )^j (8) 

2.4 Three Temperature Model for Ions, Electrons, and Radiation 

In the 3T model, the temperature in Eq. ([ 3 ]) is replaced by temperatures for ions, electrons, and radiation, 
represented by three coupled differential equations, 
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where T,T e , and T r are the ion, electron, and radiation temperatures respectively. Heat sources include 
energy deposition to the electrons due to the laser heating, Si aser , nuclear fusion particle energy deposition 
E p (negligible in these simulations), heat conduction fcVT, adiabatic compression heating and non-adiabatic 
compression heating of the ions due to artificial viscosity, q, and the viscous heating, proportional to 7y 0 Vzi. 
The electron conduction is not flux limited but the heat flux by classical conduction will be compared to the 
flux limited values as discussed in the Results. 

The deuterium-deteurium nuclear fusion reactions evolve by 
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where n B is the deuterium number density. 

Radiation transport coefficient: Rossland Opacity. The radiative diffusion coefficient, n r , in Eq. © 
is based on an approximation to the grey-body diffusive flux, T r , for the radiation energy density, E r = 
(4cr SB /c)T r 4 , as 
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where the Rossland transport mean free path is Ar ~ 1 /(pun) and kr is the Rossland averaged opacity. 


The Rossland opacity is based on a simple fit approximating the tabular opacity data 32 


2.5 Temperature Coupling Coefficients 


Ion and electron energy is exchanged through ion-electron collisions 37 38 . Electron-radiation coupling 


terms represent energy exchange between the plasma and radiation by electron-photon Compton scattering 
and by a Planck opacity absorption and emission. The coupling coefficients used here were developed in a 
simple model 32 to study a DT run-away burn problem in an infinite space 33 . It was found this simple 


model agrees with analysis and with more sophisticated grey-body diffusion computational models 32 


Ion-Electron Coupling. The ion-electron coupling term characterizes Coulomb collisions and assumes ion 
and electron distributions are near Maxwellian 37 . The ion-electron energy exchange can be characterized 


by an energy equilibration time for each ion species 38 
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where eo is the permittivity of free space, nii and m e are the masses of the ions and electrons respectively, qi 
and q e are the charges of the ions and electrons, and L is the Coulomb logarithm. The ion-electron coupling 
coefficient Ui e is then given by 


Wie = — (15) 
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where rq is the ion number density, and the ion contributions are summed over ion species. 

Radiation Coupling-Compton Scattering. The Compton scatter coupling term w c represents the transfer 
of energy from radiation to electrons through photon on free electron collisions. The coupling term is given 
in a grey-body approximation by: 
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where ctr is the Thomson cross-section, osb is the Stefan-Boltzann constant, m e is the mass of the electron, 
n e is the atomic density of electrons, and c is the speed of light. 

Radiation Coupling-Planck Opacity. The Planck opacity coupling coefficient lo p determines the rate 
at which radiation couples to the electrons by emission and absorption processes. The linearized Planck 
coupling coefficient is calculated as: 


Wp — ^pKpCTSBiTe + T 2 )(T e + T r ). 


(17) 


where the Planck opacity, K p , is evaluated as a power law fit 32 to tabulated values from the SESAME 
data base at LANL and applied here to the deuterium fuel. 


2.6 Viscous Stress Tensor in Spherical Coordinates 

The viscous stress tensor was originally evaluated as (rj 0 \7u)\7 ■ u as in Eq. © however this was re-evaluated 
in spherical coordinates to simplify in ID with the implied symmetries, iig = = 0 and d/dO = d/d(f> = 0. 

The radial stress as a function of the symmetric strain components becomes 

(V : r) r [l D] = -4^ + Tee+T ^ (18) 

r z or r 

Simplifying the stress components by the symmetry assumptions, then T rr = —2r/ 0 (du/dr) + (2/3)rj 0 V ■ u 
and Tgg = T 00 = —2r) 0 (u/r) + 2/3?y 0 V • it, with it, the radial velocity component, the tensor radial component 
becomes, 


5 






















(19) 


(V : r) r [lD] = 


+ d_ 
r 2 dr 



V • u 
3 


4 rj 0 f u V ■ it\ 
r \r 3 ) 


These four terms are interpreted as a diffusion of velocity, a compressible correction to the diffusion, a spher¬ 
ical metrics correction and the compressible contribution to the spherical coordinate correction. Evaluating 
the divergence in spherical coordinates, V-u = (1 /r 2 )d(r 2 u)/dr = du/dr+2u/r, the terms can be rearranged 
to simplify as 
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The form for viscous dissipation in the energy equation with the assumed radial symmetry, [ug = = 

d/d6 = d/d(f> = 0 ], becomes 
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Numerical issues associated with the discretization of the viscous terms are discussed briefly in the next 
section. 


3 Numerical Methodology 


The numerical scheme is based on a Lagrange velocity evaluated at mesh cell faces or nodes, i, located between 
thermodynamic quantity cell centers indexed as i and i — 1. At initial times a grid with A r sa 2 — 3/xm allows 
the Lagrange cell sizes to remain manageable (A r m i n ss 0.3/im) at the maximum compressions in this study. 
With the pressures and artificial viscosity [l7], < 7 ,; , evaluated at cells, i, and i — 1, a new (advanced time) 
Lagrange velocity, u " +1 is first determined at velocity node, i. This is used to calculate the new position of 
each Lagrangian co-ordinate and advance the density, momentum and temperatures in the Lagrange frame. 
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The three temperature equations are then solved explicitly for the viscous terms and for half the PdV 
work done by pressure on each species, 
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The three temperatures are then updated implicitly for the remaining half of the PdV work, and lastly for 
the heat conduction of each species denoted as, 


Ti 


n+1 

all 


Te 


n+1 

all 


TrT 


= V • Ki 

= v • He 
= V • K r 


A r 
A r 
A r 


(31) 

(32) 

(33) 


The implicit differencing of the diffusion operators in the computations, Eq. plj [32j [33| are performed 
using the Thomas algorithm solving a tridiagonal matrix. Following Eq. (33), the three temperature coupling 
is performed through an iterative procedure until convergence for each of the three temperatures is achieved. 

After the three temperature equations are evaluated, the species mass mixing algorithm is executed. The 
code has the option of performing mass species mixing implicitly when using only the molar gradient term 
in Eq. Q or calculating mix explicitly if other terms are desired. The species mass transport coefficient, 
p, = {Pia/p)(Pj / Vij) is computed as the harmonic mean of the cell centered values adjacent to the interface 
of interest. The dr factor from the gradients in each term is also included in the harmonic mean to account 
for the time varying Lagrange geometry of the cells. 
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(34) 


It is noted that p :j j v.j :j results in the cancellation of the density of the heavier species, n 3 , and thus avoids 
dividing by zero, where the heavy ion molar density, n 3 is zero. 

The final step of the hydrodynamics at a new time is the computation of the pressures for the ions, elec¬ 
trons, and radiation, P,. P ei P r . These are given as functions of the new values for mix species concentrations 
and temperature at each cell, 
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where tidd and tich refer to the ion number densities for the deuterium and for the CH. The electron density 
is a function of the ionization and ion density at the specified grid point, n e = (zdotidd) + ( zcrincH ) where 
z is the ionization, and is set to the atomic charge when the first shock passes. 


Viscous Terms. A centered difference approximation for the first term in Eq. (19) leads to the standard 
diffusion discretization, which is implicitly stable or can be explicitly stabilized with a time step restriction 
of order, At < A r 2 /D, where D « 77 / p. Numerical stability occurs when the discretization of the viscosity in 
the momentum equation leads to a negative coefficient for the radial velocity at the centered node, i. which in 
turn leads to a positive contribution to the diagonal coefficient in the implicit matrix solver for the viscosity. 
This stabilization occurs naturally for the fourth term in the viscous tensor, — (4/z/3r)(u/r). The second and 
third terms could be stabilized by a type of ’upwind’ differencing, where the difference approximation for 
the du = A u is written as Art ss u[i\ — u[i — 1] or Ait « u[i + 1] — u[i\, chosen to retain a negative coefficient 
for u[i\. However, this introduces numerical dissipation, and we have found the second and third terms can 
be center-differenced for higher order accuracy and are well behaved in an implicit solution method. When 
solved implicitly these new terms modify the off-diagonal coefficients without modifying the diagonal. 

For the cell centered viscous energy dissipation in the standard Lagrangian staggered grid, the velocities 
and radii values at cell faces allow a straight forward difference approximation to the velocity gradients 
evaluated in the cells. The term, (u/r), is evaluated as an average of the cell’s face values. 
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4 Results 


Basic Tests. The hydrodynamic and plasma physics models underwent extensive testing first as a one 
temperature code [34] and then in the 3T version. The code recovers the expected solutions for several 
standard problems (e.g., Sod shock tube, Guderley spherical test, planar diffusion). Grid convergence tests 
indicate as expected that the inviscid solutions do not converge in the zones at the center of symmetry for 
shock solutions imploding in spherical geometry. However, the dissipation introduced by the plasma transport 
terms allows a converged solution for the central peak temperature for grids with at least 200 zones in the 
IGF problem 34 . The plasma thermal conduction is the most important term for grid convergence of the 


temperature but plasma viscosity also plays a signficant role. The 3T temperature model was extensively 
tested 32 against other codes and in analysis for a test problem representing a DT runaway burn problem 


in an infinite space 33 


The base case. An ICF implosion of a plastic, CH, capsule of 25 micron thickness with inner radius of 400 
microns imploding on a DD fuel was selected and used to validate the ID code model and algorithms. The 
laser power is 22 TW with an absorption fraction (default set to 60%) delivered to the electrons in the outer 
CH zones over a 1 ns square pulse assumed to be representative of an Omega-scale ICF implosion. Setup 
details were previously described 34 and ICF test results from the code were found to be in reasonable 
agreement with the HELIOS ICF code 35 39 40 . The binary mix model assumes average ions in each 


material, with zch = 3.5 and Aqh = 6.5 in the CH plastic capsule. 


4.1 Comparison between Inviscid and Viscous Plasma Transport 

A radius-time plot of the Lagrangian fuel and plastic zones for the base case is shown in Fig. [ljfor the inviscid 
(a) and for the viscous (b) solutions. Each line represents a Lagrangian coordinate showing the movement 
and compression of the fuel and capsule. Laser heating of the plastic zones can be seen at the beginning of 
the simulation and the compression of the fuel continues until it is approximately 55 - 60 microns in radius. 
The shock first arrives on center between 1.4 to 1.5 ns and the return shock arrives at the inner capsule 
surface at approximately 1.7 ns, while the time of maximum compression is ss 1.9 ns. 

The incoming shock front is sharp in the inviscid case and remains steep after convergence on center. The 
first and second reflections at the fuel-capsule interface are apparent in the figure. The viscous shock has a 
less distinct front location as it converges on center and after reflection it becomes even more indistinct. The 
difference in the timing of incidence of the main shock is approximately 0.1 ns. The inviscid shock arrives 
at 1.4 ns and returns to the center point at 1.8 ns. The viscous shock arrives at about 1.5 ns with a leading 
’toe’ blurring the main shock front by about 0.1 ns. This shock also returns to the center point at about 1.8 
ns, when it appears to be much more diffuse than the inviscid case. 

As seen in Fig. [T] the fuel reaches a minimum radius of approximately 50 - 60 /irn for the inviscid and 
viscious cases The maximum compression is less for the viscous case, although the difference appears to be 
small in the radius-time plots. This will be examined in more detail when comparing the position of the 
composition fraction profiles for the fuel-capsule mixing layer. 

As the first shock is reflected from the center and returns to the capsule inner shell the inviscid shock 
is clearly sharper in the fuel and this creates a steeper shock rebounding into the plastic compared to the 
viscous case. However, Figure [I] suggests that viscous effects do not greatly affect the plastic. As the laser 
deposited energy propagates into the fuel, the plastic temperature decreases, and the strong temperature 
dependence of the ion viscosity reduces its role in the plastic. The higher average atomic mass of the plastic 
and especially the greater average ionization state also reduces the viscosity in the plastic. 
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Figure 1: (a) Radius-Time Plot for ICF Base Case without viscosity, (b) Radius-Time plot of the same case 
with viscosity. The Red and Green colors refer to fuel and plastic capsule, respectively. 



























































































































A snapshot of three temperature profiles at 1.0 ns, are shown in Fig. [2] for three cases: viscous, inviscid 
and momentum viscosity only (with no viscous dissipation in the energy equation). At this early time the 
shock fronts are near a radius of 200/zm, about half the initial fuel radius, so that convergence effects are 
still minimal. The electrons and radiation temperatures are nearly identical in each case while the ions are 
shocked to a higher temperature with a steep front in each case, and there is evidence of a small ’toe’ forming 
at the shock front in the viscous case. Electron thermal conductivity diffuses the electron temperature behind 
the shock front, and the temperatures are still small enough at this time that the electron temperature has 
not propagated beyond the ion shock front. Since the thermal conductivity for the ions is smaller by a factor 
of a Jme/rrii, the ion shock front lacks the diffusive characteristic. 

The viscous shock leads the inviscid shock at this time by about 15 cm or about 7% of the radius while 
the shock in the third case, with momentum viscosity only and no viscous dissipation in the energy equation, 
lags the inviscid case by about the same amount. The differences in shock positions are attributed primarily 
to shock energy differences with the heat from the viscous dissipation increasing the energy in the shock. 
With no viscous heating in the temperature equation, the viscosity in the momentum equation dissipates 
the shock velocity slowing the shock front compared to the inviscid case. 

A point of interest arises in comparing the shock fronts at 1 ns, in Fig. 2 to the shock fronts and arrival 
time at the target center as seen in the radius-time plots in Fig. 1. Clearly, the viscous shock leads the 
inviscid shock by ss 15 microns or 7% of the radius at 1 ns. However, once the shocks converge on center 
as seen in Fig. 1, the shock with viscosity has blurred and the main front apparently reaches the center at 
about 1.5 ns, lagging the inviscid case arriving on center at 1.4 ns. This suggests that the viscous effects vary 
significantly with the radial convergence. When there is little convergence at 1 ns the viscous shock leads the 
inviscid while the shock dissipation by viscous effects is enhanced in the convergent geometry, especially near 
the center point, leading to the viscous shock front lagging the inviscid shock during convergence at the center. 



Position (n m) 


Figure 2: Electron, Ion, and Radiation Temperature profiles for 3 cases: Inviscid, Viscous, and viscosity in 
momentum equation only labelled ’Only Viscous M’ (no viscous energy dissipation) at t — 1.0 ns. 
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A portion of the radial profiles of the temperatures are shown in Fig. [3ja) during shock convergence 
on the target center for the viscous and the invisicid cases, respectively at times of 1.52 ns and 1.42 ns. 
The inviscid and viscous profiles are compared at a time shortly after convergence (1.56 ns for each case) 
in Fig. §b). Considering first the time of shock convergence in Fig. [3ja) , the inviscid case shows a strong 
increase in temperatures at the radial center driven by the shock heating of the ions. At this time the viscous 
case shows no significant local increase in temperatures apparently due to the viscous dissipation. The ion- 
electron and electron-radiation temperature differences are significant and the differences are comparable 
between the viscid or inviscid cases. Immediately after the shock convergence in Fig. [3jb) , the central 
temperatures relax in the inviscid case to the values seen in the viscous case. 




(a) (b) 


Figure 3: Ion, Electron, and Radiation Temperatures at times (a) of first shock convergence and (b) just 
after first shock convergence. 
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Position ( ij , m) 

Figure 4: Ion, Electron, and Radiation Temperatures at times near maximum compression and neutron 
production. 


Temperature profiles in Fig. [4j are shown at a later time near maximum compression, (t = 1.85ns for 
both the inviscid and viscous cases, the time of maximum neutron production). These profiles are similar 
in shape to that seen immediately after the shock converged, although compressed in radius and slightly 
elevated in maximum temperature values on center. The viscous and inviscid profiles are very similar sug¬ 
gesting the burn dynamics will be similar as discussed next. At each of the times, in Fig. [3] and Fig. [4] there 
are large temperature differences between T), T e , T r , suggesting the 3T model is important for the inviscid 
or viscous case. The radiation energy density is small but the ion and electron temperatures are influenced 
by the radiation through the temperature coupling terms described previously. The significant temperature 
differences at these times are in contrast to the temperatures at the early implosion time, t = Ins seen in 
Fig-! when the electron and radiation temperatures are tightly coupled, and only the ion temperature is 
significantly greater due to the incoming shock heating. 


The calculation of burn weighted ion temperatures, illustrated in Fig. [5j was used as a diagnostic to 
compare the code to experimental and other numerical results. For each case, the time dependent burn 
weighted ion temperatures exhibit two local maximum values, corresponding to the times of shock heating 
and compression burn. The maximum burn weighted ion temperature in the inviscid case was found to be 
approximately 7 keV (and briefly approaches 10 keV at first shock heating), while the viscous case maximum 
burn weighted ion temperature is about 3 keV, at first shock or at maximum compression. The difference 
in time dependent burn weighted temperature between viscid and inviscid cases is large only during the 
interval, 1.44 — 1.58ns, as the first shock converges and reflects from the center. The time averaged burn 
weighted temperature in either case was similar, about 2.6 keV. These values are comparable to experimental 
results and other ID ICF simulations 


40 43 44 


Neutron production rates verses time in Fig. [6]show yield rates and yield of order 10 
10 11 and are comparable to experiment and previous simulations 


3 x 

45H47I. For each case, neutron production 
rates exhibit a local maximum value at the time of compression burn. During shock heating, neutron 
production is about 10% of the maximum level for the inviscid case and significantly less for the viscous 
case. Maximum neutron production rates were found to occur approximately at 1.85 ns, which agrees well 


with previous studies ( 24 39 40 43 , see also 45 and references therein). 


The plots of burn weighted temperatures in Fig. [5] and the neutron production rate in Fig. [6] illustrate 
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the inviscid and viscous differences in shock timings. The inviscid case burn temperature has a sharp peak 
at the time of incidence of the first shock while the viscous case does not since it dissipates the shock energy 
more effectively. The decrease in temperature indicated in Fig. [5] for the viscous case leads to a delayed 
profile for the neutron production rate compared to the inviscid simulation. The relatively small neutron 
production rate at first shock, even for the inviscid case is attributed to the lower density at that time before 
maximum compression. 



Figure 5: Burn weighted ion temperature over time for both cases, viscid and inviscid. The temperature is 
consistently higher in the inviscid case. 
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Figure 6: Neutron production rate over time for both simulations, viscid and inviscid. The peak occurs at 
the time of maximum compression as indicated in the previous plots. 


To better match experimental results during the early stage of ICF implosions it has been proposed in 
previous work to use non-local thermal transport models or to use flux-limited electron transport l](5}{7>36 . 
Figure [7] shows a plot of the classical transport electron heat flux, « e VT e in our simulations, along with the 
heat flux at the maximum physical flux limit and with the commonly implemented flux limiter of 7% of the 
maximum flux limit 16 [|T] . The maximum limited heat flux is the magnitude of the internal energy that can 
be moved by a characteristic velocity, v = v t h■ Each is plotted at a time, t = Ins, when the laser drive 
has just turned off and an equal amount of laser energy has been deposited in both the viscid and inviscid 
simulations. 
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(a) Inviscid (b) Viscous 

Figure 7: Plot of the heat flux, K e VT e , denoted (K grad T), with the kinetic flux limit, 3nkf,T e v, and the 
commonly implemented flux limit, 7% of 3 nkf,T e v 


The comparison for the viscous and inviscid cases illustrates a few points. The shape of the heat flux 
front is considerably steeper in the inviscid case, and the difference in the size of the shocks is evident. In 
the inviscid case, the classic heat flux exceeds the 7% limited flux over the radius, 195 — 210 /zm while in 
the viscous case, this occurs over the radius 180 — 205 These are the regions which might require the 
implementation of a flux limiter, and the region is « 60% larger in the viscous case. This suggests that 
viscous effects can modify the role of a flux limiter. Where these flux limiters are implemented with phe¬ 
nomenological model parameters, the parameter values may need adjustment when including plasma viscous 
effects. Conversely, including viscosity in simulations may reduce the need for phenomenological adjustment 
of the electron heat flux. 


Figure [8] shows the mass composition profiles due to plasma species transport across the fuel-capsule 
interface at the time of maximum compression, t ~ 1.85ns. The radial locations of the midpoints of the 
mixing layers at this time are 55 microns for the inviscid case, and 58.5 microns for the viscous case, where 
the difference of approximately 3.5 /im is 6.4% of the inviscid radius and corresponds to a difference of 20% in 
fuel volume. The mixing layer itself is slightly wider in the viscous case, ~ 2.5/xm, than in the inviscid case 
at ss 2/i?n. With the differences in mix width and radius taken together this corresponds to a mixed region 
volume about 40% larger in the viscous case. A ID HELIOS simulation of similar ICF implosions 24 showed 
greater convergence (r m i„ ~ 20 but post-processed calculations from those simulations showed good 
agreement with the mass transport mix widths in this study, suggesting comparable density and temperatures 
near the fuel-capsule interface in the two cases. The convergence differences can be related to laser input 
and initial fill pressure differences. The widths of these mix layers by plasma gradient driven transport is 
small in these idealized ID simulations, however, we expect the mass diffusivity will play an enhanced role in 
2D and 3D 25j|26 where the fuel-capsule interfacial area can be increased significantly by instability growth 
and possibly by turbulence. 
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t = 1.854 



x m) x (/i m) 


(a) Inviscid 


(b) Viscous 


Figure 8: Plot of the mixing layer between DD and CH at the time of compression burn for (a) the Inviscid 
Simulation and (b) Viscous Simulation 


We considered results in Ref. [3] as a guide for expected asymmetry of the mass species mixing profiles. 
Our profile results are more symmetric between the fuel and capsule than expected, suggesting that the 
carbon-based capsule may not be sufficiently high Z, to achieve the asymmetry predicted in 13 . 

Results Summary. Select results from the inviscid and viscous simulations to 2.5 ns are summarized 
in Table I. The values indicate for the inviscid case that both maximum temperatures and densities at the 
origin occur during the ’spike’ at first shock convergence (t ~ 1.43ns) leading to high pressures in that 
spike. The viscous case maximum densities and temperatures at the origin occur during compression and at 
slightly different times, 1.83 ns and 2.0 ns respectively. The offset in time between the maxima in density 
and temperature leads to a greatly reduced instantaneous pressure maximum at the origin. 


Table 1: ICF Result Summary 


Model Parameter 

Inviscid case 

Viscous case 

Final time -ns- 

2.50 

2.50 

Highest temperature (ion) at origin 

10392.35 

4177.48 

Highest temperature (e) at origin 

1689.5 

1579.69 

Highest temperature occurs at time -ns- 

1.430 

1.833 

Highest density at origin 10 1A cm~ 6 

128.95 

3.216 

Highest density occurs at time -ns- 

1.430 

2.027 

Maximum pressure at origin (Gbar) 

117.31 

1.4139 

Ion Temp in Zone 1 at final time (eV) 

405.06 

433.72 

Elec Temp in Zone 1 at final time (eV) 

415.67 

446.37 

Rad Temp in Zone 1 at final time (eV) 

425.91 

454.92 

Ion Temp in Outer Fuel Zone at final time (eV) 

152.7 

162.2 

Elec Temp in Outer Fuel Zone at final time (eV) 

152.6 

162.1 

Rad Temp in Outer Fuel Zone at final time (eV) 

152.7 

162.3 

Burn Weighted Ion Temperature (eV) 

2704. 

2655. 


At the final run time, 2.5 ns, the three temperatures at the origin are similar but not completely equili¬ 
brated in the inviscid or viscous case. The temperatures are slightly greater in the viscous cases, as expected 
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due to the viscous heating. At this late time, in the outer fuel zone at the fuel-plastic boundary, the three 
temperatures have equilibrated in the inviscid and viscous cases with a slightly higher final temperature 
for the viscid case. The time averaged burn weighted ion temperatures are nearly identical for the viscid 
and inviscid cases, in sharp contrast to the maximum temperature values which spiked at the origin for the 
inviscid case. 

Artificial viscosity. Trade-offs between physical plasma viscosity and artificial viscosity in Lagrange 
numerical methods have been examined in 27 . Our results compliment and extend their findings. We 


found that with plasma viscosity in the simulation, the solution remains numerically stable if the artificial 
viscosity is turned off after the shock converges on the target center. After this time, the entire domain has 
become a fully ionized plasma with significant viscosity. In the simulations including plasma viscosity and 
with artificial viscosity turned off, the solutions show considerably more fluctuations compared to the inviscid 
or the viscous case when artificial viscosity is included. These fluctuations are seen in the fuel and capsule 
regions especially as the shocks pass from the center back into the capsule material. It is hypothesized that 
the increased fluctuations, apparent in the absence of artificial viscosity, lead to increases in kinetic energy 
and entropy which can reduce implosion efficiency. This result was observed near the end of the current 
project and is planned to be explored in future work. 


5 Conclusion 


In this study, a one dimensional, staggered grid Lagrangian hydrodynamic code was written to include three 
temperature effects for ions, electrons, and radiation with plasma transport and dissipation including viscous 
diffusion of momentum, viscous heating of ions, and species mass transport, in addition to the usual plasma 
physics describing species thermal conductivities and temperature coupling through collisional equilibration. 
This code was used to examine ICF implosions typical of certain DD-CH fuel-capsules at the Omega Laser 
Facility. 

The goal was to show first that the code recovers reasonable and realistic results compared to experiments 
and other codes for this class of ICF implosions while then examining the differences in result details for 
A-B comparison simulations which differ only in the plasma viscosity. In each of these cases, the incidence 
of the first shock on target center occurs at 1.4 - 1.5 ns in the simulations and the return shock arrives at 1.8 
- 1.9 ns. These shock timings are in generally good agreement with previous calculations and measurements 
IIIHIL 43 ]. The neutron production rate, and the maximum and burn weighted ion temperatures 
match expected results 44-46 . The minimum radius of the fuel during compression reached « 50/xm, in 
contrast to previous calculations of ss 20 fj,m using HELIOS 24,39 40 , which can be related to differences 


in the prescribed laser inputs, thickness of the capsule, and fuel fill pressures. 

Our examination of viscosity in ICF implosions, shows the plasma viscosity and viscous heating of the 
ions resulted in differences in maximum fuel compression, neutron production rate, and the time dependent 
burn weighted ion temperatures. The viscous shock front is about 8% ahead of the inviscid shock at the 
early time (1 ns) position of the incoming shocks. While the inviscid shock is slower at 1 ns, it arrives earlier 
on axis by approximately 0.1 ns. This behavior may be related to the radial convergence terms in the viscous 
dissipation which play a minor role at the large radius at early time, but play a more important role as the 
convergence becomes significant at the times of first shock collapse on target center and during maximum 
compression. 

The fuel compresses to a minimum radius for the viscous case that is greater by 3.5 /rm, or 6%, compared 
to the inviscid compression. This corresponds to a 20% increase in fuel volume for the viscous case, and 
combined with a slightly wider mix width for the viscous case there is a 40% increase in the mix volume 
for the viscous case. The widths for the mass composition mixing layer of a few microns calculated from 
post-processed HELIOS simulations 24j391 are in general agreement with the viscous and the inviscid results 
here. It is evident from the results that viscous effects cause a decrease in the compression. This increase in 
the size of the compressed region with viscosity reduces the rho-r of the fuel and is also related to an increase 
in entropy deposited in the fuel and capsule due to the viscous dissipation through both the momentum and 
energy equations. It is well known that entropy production degrades confinement and is a key concern in 
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implosion efficiency 36 


The ion temperatures are consistently higher for the inviscid case at first shock convergence. This leads 
to a factor of approximately two difference in the ion temperatures at that time and it follows that the 
neutron production rate is also augmented during shock convergence. The time averaged burn weighted ion 
temperatures are only slightly larger in the inviscid case and at maximum compression they agree within 
a few percent. The neutron yields are comparable during this phase between inviscid and viscous cases. It 
was also determined that the implementation of viscous effects modifies the need for an electron heat flux 
limiter in these hydrodynamic simulations. As the viscous terms increase the width of the shock front and 
decrease the magnitude of the ion temperature in the shocked region, the need for a flux limiter may be 
reduced across the domain with viscosity included. 

After the first shock converges, the entire region becomes a plasma with significant temperatures driving 
up the magnitude of the plasma viscosity. After that time, the plasma viscosity and viscous dissipation 
eliminate the need for artificial viscosity to stabilize the numerics, and the artificial viscosity can be turned 
off. Results suggest that the subsequent solutions allow greater shock-wave fluctuations in the fuel and 
especially in the capsule material compared to cases when artificial viscosity is included. Previous work in 
this area 


27 confirms this is a topic worth examining further and in greater detail. 


It was determined that viscous effects play a crucial role in implosion dynamics for this DD-CH model 
suggesting that other codes should consistently include viscous effects to model ICF implosions with greater 
fidelity. Viscosity influences shock timings, fuel compression, neutron production rates, and heat fluxes 
especially at the time of first shock convergence. Viscous effects during the later compression phase appear 
to be small in this ID example, but are likely to play a more important role in 2D and 3D as suggested 
in previous studies 25 26 . It is expected that viscosity and plasma species mass transport will have a 


significant influence in other ICF capsule implosions. 
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